/*
 * com.javable.cese -- Cell Electrophysiology Simulation Environment (CESE)
 *
 * Copyright (C) 1999 - 2005 Sergey Missan
 *
 * Dalhousie Univ., Halifax, NS, Canada
 *
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public
 * License along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
 *
 */

package models.com.javable.cese.models;

import java.io.Serializable;

import com.javable.cese.templates.AbstractAPModel;

/**
 * The Luo-Rudy I Model of Mammalian Ventricular Cardiac Action Potentials,
 * 1991. Detailed list of equations and model description are provided in
 * <p>
 * Circ Res 1991;68:1501-1526.
 */

public class LuoRudyModel_I_1991 extends AbstractAPModel implements
        Serializable {

    /*
     * The beginning here lists all the variables and defines them (this code
     * uses all global variables)
     */

    /* Cell Geometry */
    private double l = 0.01; // Length of the cell (cm)

    public synchronized void setl(double sl) {
        l = sl;
        getCellGeometry();
    }

    public double getl() {
        return l;
    }

    private double a = 0.0011; // Radius of the cell (cm)

    public synchronized void seta(double sa) {
        a = sa;
        getCellGeometry();
    }

    public double geta() {
        return a;
    }

    private double vcell; // Cell volume (uL)

    public double getvcell() {
        return vcell;
    }

    private double ageo; // Geometric membrane area (cm^2)

    public double getageo() {
        return ageo;
    }

    private double acap; // Capacitive membrane area (cm^2)

    public double getacap() {
        return acap;
    }

    private double vmyo; // Myoplasm volume (uL)

    public double getvmyo() {
        return vmyo;
    }

    private double vmito; // Mitochondria volume (uL)

    public double getvmito() {
        return vmito;
    }

    private double vcleft; // Cleft volume (uL)

    public double getvcleft() {
        return vcleft;
    }

    /* Terms for Solution of Conductance and Reversal Potential */
    public static final double R = 8314; // Universal Gas Constant (J/kmol*K)
    public static final double frdy = 96485; // Faraday's Constant (C/mol)

    private double temp = 310; // Temperature (K)

    public synchronized void settemp(double stemp) {
        temp = stemp;
    }

    public double gettemp() {
        return temp;
    }

    /* Ion Valences */
    public static final double zna = 1; // Na valence
    public static final double zk = 1; // K valence
    public static final double zca = 2; // Ca valence

    /* Ion Concentrations */
    private double nai = 18.0; // Intracellular Na Concentration (mM)

    public synchronized void setnai(double snai) {
        nai = snai;
    }

    public double getnai() {
        return nai;
    }

    private double nao = 140.0; // Extracellular Na Concentration (mM)

    public synchronized void setnao(double snao) {
        nao = snao;
    }

    public double getnao() {
        return nao;
    }

    private double ki = 145.0; // Intracellular K Concentration (mM)

    public synchronized void setki(double ski) {
        ki = ski;
    }

    public double getki() {
        return ki;
    }

    private double ko = 5.4; // Extracellular K Concentration (mM)

    public synchronized void setko(double sko) {
        ko = sko;
    }

    public double getko() {
        return ko;
    }

    private double cai = 0.0002; // Intracellular Ca Concentration (mM)

    public synchronized void setcai(double scai) {
        cai = scai;
    }

    public double getcai() {
        return cai;
    }

    private double cao = 1.8; // Extracellular Ca Concentration (mM)

    public synchronized void setcao(double scao) {
        cao = scao;
    }

    public double getcao() {
        return cao;
    }

    /* Myoplasmic Na Ion Concentration Changes */
    private double naiont; // Total Na Ion Flow (mM/ms)

    public double getnaiont() {
        return naiont;
    }

    private double dnai; // Change in Intracellular Na Concentration (mM)

    public double getdnai() {
        return dnai;
    }

    /* Myoplasmic K Ion Concentration Changes */
    private double kiont; // Total K Ion Flow (mM/ms)

    public double getkiont() {
        return kiont;
    }

    private double dki; // Change in Intracellular K Concentration (mM)

    public double getdki() {
        return dki;
    }

    /* Myoplasmic Ca Ion Concentration Changes */
    private double caiont; // Total Ca Ion Flow (mM/ms)

    public double getcaiont() {
        return caiont;
    }

    private double dcai; // Change in myoplasmic Ca concentration (mM)

    public double getdcai() {
        return dcai;
    }

    /* Fast Sodium Current (time dependent) */
    private double ina; // Fast Na Current (uA/uF)

    public synchronized void setina(double sina) {
        ina = sina;
    }

    public double getina() {
        return ina;
    }

    private double gna; // Max. Conductance of the Na Channel (mS/uF)

    public synchronized void setgna(double sgna) {
        gna = sgna;
    }

    public double getgna() {
        return gna;
    }

    private double ena; // Reversal Potential of Na (mV)

    public synchronized void setena(double sena) {
        ena = sena;
    }

    public double getena() {
        return ena;
    }

    private double ah; // Na alpha-h rate constant (ms^-1)

    public double getah() {
        return ah;
    }

    private double bh; // Na beta-h rate constant (ms^-1)

    public double getbh() {
        return bh;
    }

    private double aj; // Na alpha-j rate constant (ms^-1)

    public double getaj() {
        return aj;
    }

    private double bj; // Na beta-j rate constant (ms^-1)

    public double getbj() {
        return bj;
    }

    private double am; // Na alpha-m rate constant (ms^-1)

    public double getam() {
        return am;
    }

    private double bm; // Na beta-m rate constant (ms^-1)

    public double gebm() {
        return bm;
    }

    private double h; // Na inactivation h

    public double geth() {
        return h;
    }

    private double j; // Na inactivation j

    public double getj() {
        return j;
    }

    private double m; // Na activation m

    public double getm() {
        return m;
    }

    private double scgna = 1;

    // Scaling factor for Max. Conductance of the Na Channel
    public synchronized void setscgna(double sscgna) {
        scgna = sscgna;
    }

    public double getscgna() {
        return scgna;
    }

    /* Inward Ca Current */
    private double ica; // Total current through the Inward Ca channel (uA/uF)

    public synchronized void setica(double sica) {
        ica = sica;
    }

    public double getica() {
        return ica;
    }

    private double gca; // Max. Ca current through Ca channel (uA/uF)

    public synchronized void setgca(double sgca) {
        gca = sgca;
    }

    public double getgca() {
        return gca;
    }

    private double eca; // Reversal Potential of Inward Ca Current (mV)

    public synchronized void seteca(double seca) {
        eca = seca;
    }

    public double geteca() {
        return eca;
    }

    private double ad; // Ca alpha-d rate constant (ms^-1)

    public double getad() {
        return ad;
    }

    private double bd; // Ca beta-d rate constant (ms^-1)

    public double getbd() {
        return bd;
    }

    private double d; // Voltage dependent activation gate d

    public double getd() {
        return d;
    }

    private double af; // Ca alpha-f rate constant (ms^-1)

    public double getaf() {
        return af;
    }

    private double bf; // Ca beta-f rate constant (ms^-1)

    public double getbf() {
        return bf;
    }

    private double f; // Voltage dependent inactivation gate f

    public double getf() {
        return f;
    }

    private double scgca = 1;

    // Scaling factor for Max. Conductance of the Inward Ca channel
    public synchronized void setscgca(double sscgca) {
        scgca = sscgca;
    }

    public double getscgca() {
        return scgca;
    }

    /* Time-dependent Potassium Current */
    private double ikt; // Time-dependent K Current (uA/uF)

    public synchronized void setikt(double sikt) {
        ikt = sikt;
    }

    public double getikt() {
        return ikt;
    }

    private double gkt;

    // Channel Conductance of Time-dependent K Current (mS/uF)
    public synchronized void setgkt(double sgkt) {
        gkt = sgkt;
    }

    public double getgkt() {
        return gkt;
    }

    private double ekt; // Reversal Potential of Time-dependent K Current (mV)

    public synchronized void setekt(double sekt) {
        ekt = sekt;
    }

    public double getekt() {
        return ekt;
    }

    private double axt; // K alpha-xt rate constant (ms^-1)

    public double getaxt() {
        return axt;
    }

    private double bxt; // K beta-xt rate constant (ms^-1)

    public double getbxt() {
        return bxt;
    }

    private double xt; // Time-dependent K activation

    public double getxt() {
        return xt;
    }

    private double r; // K time-independent inactivation r

    public double getr() {
        return r;
    }

    private double scgkt = 1;

    // Scale factor for Channel Conductance of Time-dependent K Current
    public synchronized void setscgkt(double sscgkt) {
        scgkt = sscgkt;
    }

    public double getscgkt() {
        return scgkt;
    }

    /* Potassium Current (time-independent) */
    private double iki; // Time-independent K current (uA/uF)

    public synchronized void setiki(double siki) {
        iki = siki;
    }

    public double getiki() {
        return iki;
    }

    private double gki;

    // Channel Conductance of Time Independent K Current (mS/uF)
    public synchronized void setgki(double sgki) {
        gki = sgki;
    }

    public double getgki() {
        return gki;
    }

    private double eki;

    // Reversal Potential of Time Independent K Current (mV)
    public synchronized void seteki(double seki) {
        eki = seki;
    }

    public double geteki() {
        return eki;
    }

    private double aki; // K alpha-ki rate constant (ms^-1)

    public double getaki() {
        return aki;
    }

    private double bki; // K beta-ki rate constant (ms^-1)

    public double getbki() {
        return bki;
    }

    private double kin; // K inactivation

    public double getkin() {
        return kin;
    }

    private double scgki = 1;

    // Scaling factor for Channel Conductance of Time Independent K Current
    public synchronized void setscgki(double sscgki) {
        scgki = sscgki;
    }

    public double getscgki() {
        return scgki;
    }

    /* Plateau Potassium Current */
    private double ikp; // Plateau K current (uA/uF)

    public synchronized void setikp(double sikp) {
        ikp = sikp;
    }

    public double getikp() {
        return ikp;
    }

    private double gkp; // Channel Conductance of Plateau K Current (mS/uF)

    public synchronized void setgkp(double sgkp) {
        gkp = sgkp;
    }

    public double getgkp() {
        return gkp;
    }

    private double ekp; // Reversal Potential of Plateau K Current (mV)

    public synchronized void setekp(double sekp) {
        ekp = sekp;
    }

    public double getekp() {
        return ekp;
    }

    private double kp; // K plateau factor

    public double getkp() {
        return kp;
    }

    private double scgkp = 1;

    // Scaling factor for Channel Conductance of Plateau K Current
    public synchronized void setscgkp(double sscgkp) {
        scgkp = sscgkp;
    }

    public double getscgkp() {
        return scgkp;
    }

    /* Background Current */
    private double ib; // Background current (uA/uF)

    public synchronized void setib(double sib) {
        ib = sib;
    }

    public double getib() {
        return ib;
    }

    private double gb; // Channel Conductance of Background current (uA/uF)

    public double getgb() {
        return gb;
    }

    private double scgb = 1;

    // Scaling factor for Max. conductance of Background
    public synchronized void setscgb(double sscgb) {
        scgb = sscgb;
    }

    public double getscgb() {
        return scgb;
    }

    /* ** CONSTRUCTOR ** */
    public LuoRudyModel_I_1991() {
        super();
        getCellGeometry();
    }

    public void setCurrents() {
        currentINa();
        currentICa();
        currentIKt();
        currentIKi();
        currentIKp();
        currentIb();
        currentIt();
    }

    public void setConcentrations() {
         concentrationCai();
    }

    public void doSimulate() {
        setCurrents();
        setConcentrations();

        vnew = v - it * dt;
        dvdtnew = (vnew - v) / dt;
        t = t + dt;
        v = vnew;
        dvdt = dvdtnew;
    }

    /** ***************************************************** */
    /**
     * Cell Geometry
     */
    protected void getCellGeometry() {
        vcell = 1000 * Math.PI * a * a * l;
        ageo = 2 * Math.PI * a * a + 2 * Math.PI * a * l;
        acap = ageo * 2;
        vmyo = vcell * 0.68;
        vmito = vcell * 0.26;
        vcleft = vcell * 0.12 / 0.88;
    }

    /* Functions that describe the currents begin here */

    /**
     * Calculates Fast Na Current
     */
    protected void currentINa() {
        gna = 23.0 * scgna;
        ena = ((R * temp) / frdy) * Math.log(nao / nai);

        if (v != -47.13) {
            am = 0.32 * (v + 47.13) / (1.0 - Math.exp(-0.1 * (v + 47.13)));
        } else {
            am = 3.2;
        }
        bm = 0.08 * Math.exp(-v / 11.0);

        if (v < -40.0) {
            ah = 0.135 * Math.exp((80.0 + v) / -6.8);
            bh = 3.56 * Math.exp(0.079 * v) + 310000 * Math.exp(0.35 * v);
            aj = (-127140 * Math.exp(0.2444 * v) - 0.00003474 * Math
                    .exp(-0.04391 * v))
                    * (v + 37.78) / (1.0 + Math.exp(0.311 * (v + 79.23)));
            bj = (0.1212 * Math.exp(-0.01052 * v))
                    / (1.0 + Math.exp(-0.1378 * (v + 40.14)));
        } else {
            ah = 0.0;
            bh = 1.0 / (0.13 * (1.0 + Math.exp((v + 10.66) / -11.1)));
            aj = 0.0;
            bj = (0.3 * Math.exp(-0.0000002535 * v))
                    / (1.0 + Math.exp(-0.1 * (v + 32.0)));
        }

        if (t == 0) // Steady State Values
        {
            h = ah / (ah + bh);
            j = aj / (aj + bj);
            m = am / (am + bm);
        } else {
            h = ah / (ah + bh) - ((ah / (ah + bh)) - h)
                    * Math.exp(-dt / (1 / (ah + bh)));
            j = aj / (aj + bj) - ((aj / (aj + bj)) - j)
                    * Math.exp(-dt / (1 / (aj + bj)));
            m = am / (am + bm) - ((am / (am + bm)) - m)
                    * Math.exp(-dt / (1 / (am + bm)));
        }

        ina = gna * m * m * m * h * j * (v - ena);
    }

    /**
     * Calculates Inward Ca Current
     */
    protected void currentICa() {
        gca = 0.09 * scgca;
        eca = 7.7 - 13.0287 * Math.log(cai);

        ad = 0.095 * Math.exp(-0.01 * (v - 5.0))
                / (1.0 + Math.exp(-0.072 * (v - 5.0)));
        bd = 0.07 * Math.exp(-0.017 * (v + 44.0))
                / (1.0 + Math.exp(0.05 * (v + 44.0)));

        af = 0.012 * Math.exp(-0.008 * (v + 28.0))
                / (1.0 + Math.exp(0.15 * (v + 28.0)));
        bf = 0.0065 * Math.exp(-0.02 * (v + 30.0))
                / (1.0 + Math.exp(-0.2 * (v + 30.0)));

        if (t == 0) // Steady State Values
        {
            d = ad / (ad + bd);
            f = af / (af + bf);
        } else {
            d = ad / (ad + bd) - ((ad / (ad + bd)) - d)
                    * Math.exp(-dt / (1 / (ad + bd)));
            f = af / (af + bf) - ((af / (af + bf)) - f)
                    * Math.exp(-dt / (1 / (af + bf)));
        }

        ica = gca * d * f * (v - eca);
    }

    /**
     * Calculates Time-dependent K Current
     */
    protected void currentIKt() {
        gkt = (0.282 * Math.sqrt(ko / 5.4)) * scgkt;
        ekt = ((R * temp) / frdy)
                * Math.log((ko + 0.01833 * nao) / (ki + 0.01833 * nai));

        if (v <= -100.0) {
            r = 1.0;
        } else {
            if (v != -77.0) {
                r = 2.837 * (Math.exp(0.04 * (v + 77.0)) - 1.0)
                        / ((v + 77.0) * Math.exp(0.04 * (v + 35.0)));
            } else {
                r = 0.60888;
            }
        }

        axt = 0.0005 * Math.exp(0.083 * (v + 50.0))
                / (1.0 + Math.exp(0.057 * (v + 50.0)));
        bxt = 0.0013 * Math.exp(-0.06 * (v + 20.0))
                / (1.0 + Math.exp(-0.04 * (v + 20.0)));

        if (t == 0) // Steady State Values
        {
            xt = axt / (axt + bxt);
        } else {
            xt = axt / (axt + bxt) - ((axt / (axt + bxt)) - xt)
                    * Math.exp(-dt / (1 / (axt + bxt)));
        }

        ikt = gkt * xt * r * (v - ekt);
    }

    /**
     * Calculates Time-Independent K Current
     */
    protected void currentIKi() {
        gki = (0.6047 * Math.sqrt(ko / 5.4)) * scgki;
        eki = ((R * temp) / frdy) * Math.log(ko / ki);

        aki = 1.02 / (1.0 + Math.exp(0.2385 * (v - eki - 59.215)));
        bki = (0.49124 * Math.exp(0.08032 * (v - eki + 5.476)) + Math
                .exp(0.06175 * (v - eki - 594.31)))
                / (1.0 + Math.exp(-0.5143 * (v - eki + 4.753)));

        kin = aki / (aki + bki);

        iki = gki * kin * (v - eki);
    }

    /**
     * Calculates Plateau K Current
     */
    protected void currentIKp() {
        gkp = 0.0183 * scgkp;
        ekp = eki;

        kp = 1.0 / (1.0 + Math.exp((7.488 - v) / 5.98));

        ikp = gkp * kp * (v - ekp);
    }

    /**
     * Calculates Background Current
     */
    protected void currentIb() {
        gb = 0.03921 * scgb;

        ib = gb * (v + 59.87);
    }

    /**
     * Calculates Total Current
     *
     * Total sum of currents is calculated here, if the time is between stimtime =
     * 0 and stimtime = tstim, a stimulus is applied
     */
    public void currentIt() {
        naiont = ina;
        kiont = ikt + iki + ikp + ib;
        caiont = ica;

        if ((t >= tstim) && (t < (tstim + dt))) {
            stimtime = 0;
            tstim = tstim + bcl;
            dvdtmax = false;
        }

        if ((stimtime >= 0) && (stimtime < stl)) {
            it = st + naiont + kiont + caiont;
        } else {
            it = naiont + kiont + caiont;
        }

        stimtime = stimtime + dt;
    }

    /* Functions that calculate intracellular ion concentrations begins here */

    /**
     * Calculates new myoplasmic Ca ion concentration
     */
    protected void concentrationCai() {
        dcai = dt * ((-0.0001 * caiont) + 0.07 * (0.0001 - cai));
        cai = dcai + cai;
    }

}